list_of_packages <- c("tidyverse", "plotly", "here", "signal", "DT")
new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.4
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(here)
## here() starts at /home/francisko/coding/r/flight-path-experiment
library(signal)
##
## Attaching package: 'signal'
## The following object is masked from 'package:plotly':
##
## filter
## The following object is masked from 'package:dplyr':
##
## filter
## The following objects are masked from 'package:stats':
##
## filter, poly
library(DT)
b195_t1 <- readr::read_csv(here::here("unpaired-points", "b195-t1-unpaired-points-xyz.csv"))
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## x_1 = col_double(),
## y_1 = col_double(),
## z_1 = col_double()
## )
I created the filter using the same parameter values Camille describes in his paper.
butterworth_filter <-
butter(n = 4, # order
W = 0.5, # critical frequency
type = "low")
b195_t1 %>%
tidyr::drop_na() %>%
dplyr::mutate(x_butterworth = signal::filter(butterworth_filter, x_1),
y_butterworth = signal:: filter(butterworth_filter, y_1),
z_butterworth = signal::filter(butterworth_filter, z_1)) -> b195_t1_filtered
b195_t1_filtered %>%
DT::datatable(extensions = 'Buttons',
options = list(dom = 'Blfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
lengthMenu = list(c(10,25,50,-1),
c(10,25,50,"All"))))
Here is a plot with both the raw and the “smoothed” paths so you can see difference.
plot_ly(data = b195_t1_filtered,
type = "scatter3d",
name = "Raw tracking",
x = b195_t1_filtered$x_1,
y = b195_t1_filtered$y_1,
z = b195_t1_filtered$z_1,
mode = "lines",
line = list(color = "steelblue")) %>%
add_trace(
type = "scatter3d",
name = "Tracking smoothed usind butterworth filter",
x = b195_t1_filtered$x_butterworth,
y = b195_t1_filtered$y_butterworth,
z = b195_t1_filtered$z_butterworth,
mode = "lines",
line = list(color = "firebrick")
)
I used dist(), a function from base R. It uses Euclidean distances by default and, since I did not change it, that’s the distances I sent you in that other Rmd file.
I’ll give a brief tour of how I calculated the distances to clarify things.
I defined the following function to calculate distance, it:
dist to calculate the distances;calc_flight_dist <- name <- function(dat) {
dat_dist <- {{dat}} %>%
dplyr::filter(!is.na(x_1)) %>% # removes NA
dist() %>% # calculates distance
as.matrix() # converts it to a matrix
return(sum((dat_dist[row(dat_dist) == col(dat_dist) + 1]))) # sums the distance between points by summing the points in line n and column n + 1
}
Here is it working in practice. First, I define a toy set of points.
tibble(
x_1 = c(0, 3, 3, 3), # there' only movement along this axis
y = c(0, 0, 1, 1),
z = c(0, 0, 0, 4)
) -> ex1
Here’s a quick plot so it is easier to visualize the path:
plot_ly(data = ex1,
type = "scatter3d",
name = "Raw tracking",
x = ex1$x_1,
y = ex1$y,
z = ex1$z,
mode = "lines",
line = list(color = "steelblue"))
Now let’s apply the distance calculation function:
calc_flight_dist(ex1)
## [1] 8
It returns 8, which is the total distance in this path.
To calculate sinuosity, we need the distance between the first and last points, which we can get in the distance matrix, in the last value of the first line (or the first row and last column). We get this using the function nrow (or ncol) that returns the total number of rows in the matrix and thus we can use it so get the value in the last row by using it as the index of a subset call:
dist_dat <-
b195_t1 %>%
dplyr::filter(!is.na(x_1)) %>%
dist() %>%
as.matrix()
dist_dat[1, nrow(dist_dat)] # here it is
## [1] 2.124259